%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                               Figures                                   %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% --------------------------- Description --------------------------------- 

    % This script plots some of the charts that are displayed in the 
    % Appendix and Supplementary Material of the Article

%--------------------------------------------------------------------------
%                        Supplementary Material                           %
%--------------------------------------------------------------------------

% --------------------------- Figure 2 ------------------------------------

    % Compute the expected product age given firm age and products
    ExpcAge_An = permute(sum(AgeDensity .* repmat(avec', 1, length(avec), p.maxprods), 1),[3 2 1]);

    % Sum over number of products for an average
    ExpcAge_A = sum(ExpcAge_An .* CondProbs, 1);

    % Figure 2(a)
    figure('Name','Supplementary Material - Figure 2');
    subplot(1,2,1)
    plot(avec, ExpcAge_An(1, :), 'linewidth', 2, 'DisplayName','One Product');
    hold on
    plot(avec, ExpcAge_An(2, :), 'linewidth', 2, 'DisplayName','Two Products');
    hold on
    plot(avec, ExpcAge_An(3, :), 'linewidth', 2, 'DisplayName','Three Products');
    hold on
    plot(avec, ExpcAge_An(4, :), 'linewidth', 2, 'DisplayName','Four Products');
    hold on
    plot(avec, ExpcAge_An(5, :), 'linewidth', 2, 'DisplayName','Five Products');
    hold on
    plot(avec, ExpcAge_A, '--','linewidth', 2, 'DisplayName','Average Across Firms');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Firm Age", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Average Product Age}'}, 'FontWeight', 'normal')
    title("Average Product Age by Firm Age", "FontSize", 14)
    legend('Location', 'southeast')

    % Figure 2(b)
    subplot(1,2,2)
    plot(avec, CondProbs(1, :), 'linewidth', 2, 'DisplayName','One Product');
    hold on
    plot(avec, CondProbs(2, :), 'linewidth', 2, 'DisplayName','Two Products');
    hold on
    plot(avec, CondProbs(3, :), 'linewidth', 2, 'DisplayName','Three Products');
    hold on
    plot(avec, CondProbs(4, :), 'linewidth', 2, 'DisplayName','Four Products');
    hold on
    plot(avec, CondProbs(5, :), 'linewidth', 2, 'DisplayName','Five Products');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Firm Age", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Cond. Probability Number of Products}'}, 'FontWeight', 'normal')
    title("Probability of Product Portfolios by Firm Age", "FontSize", 14)
    legend('Location', 'northeast')

%--------------------------------------------------------------------------
%                            Appendix                                     %
%--------------------------------------------------------------------------

% --------------------------- Figure 5 ------------------------------------

    figure('Name','Appendix - Figure 5');
    plot(avec, FirmAge_dens(avec), 'linewidth', 2, 'DisplayName','2% Population Growth');
    hold on
    plot(avec, FirmAge_dens_counter(avec), 'linewidth', 2, 'DisplayName','0.24% Population Growth');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Firm Age", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Density}'}, 'FontWeight', 'normal')
    title("Firm Age Distribution", "FontSize", 14)
    legend('Location', 'northeast')

% --------------------------- Figure 4 ------------------------------------

    figure('Name','Appendix - Figure 4');
    subplot(1,2,1)
    plot(avec, Exit_age*100, 'linewidth', 2, 'DisplayName','2% Population Growth');
    hold on
    plot(avec, Exit_age_counter*100, 'linewidth', 2, 'DisplayName','0.24% Population Growth');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Firm Age", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Lifecycle Exit Rate (%)}'}, 'FontWeight', 'normal')
    title("Exit Rates by Age", "FontSize", 14)
    legend('Location', 'northeast')

    subplot(1,2,2)
    plot(avec, Esales, 'linewidth', 2, 'DisplayName','2% Population Growth');
    hold on
    plot(avec, Esales_counter, 'linewidth', 2, 'DisplayName','0.24% Population Growth');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Firm Age", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Lifecycle Sales Growth (Log Points)}'}, 'FontWeight', 'normal')
    title("Sales by Age", "FontSize", 14)
    legend('Location', 'northeast')

 
% --------------------------- Figure 2 ------------------------------------

    % Determine limits for the transformed productivity grid
    a = 850;
    b = 1050;

    % Determine (end) limit for the gap grid
    c = 20;

    % Determine the grid for transformed efficiency to be used
    qgrid = p.qgrid_adj;

    % Obtain the density associated with varieties with competitors
    delta_q = qgrid(2) - qgrid(1);
    FdensC = zeros(length(Deltagrid), p.n3);
    FdensC(:, 1:end-1) = diff(F1/F1(end,end), 1, 2)/delta_q;

    % For every possible value of productivity...
    for q = 1:p.n3

        % Compute the derivative with respect to Delta
        FdensC(:, q) = gradient(FdensC(:, q), Deltagrid);

    end
    
    % Prepare data for surface plot
    [X,Y] = meshgrid(qgrid(a:b),Deltagrid(1:c));
    Z = FdensC(1:c,a:b);
    
    figure('Name','Appendix - Figure 2');
    subplot(1,2,1)
    surf(X,Y,Z)
    shading interp
    view(2)
    xlim([qgrid(a) qgrid(b)])
    ylim([Deltagrid(1) Deltagrid(c)])
    xlabel("Transformed Productivity, q_{hat}", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Gap}'}, 'FontWeight', 'normal')
    title("Joint Distribution of Efficiency and Gaps", "FontSize", 14)

   % Determine limits for the transformed productivity grid
    a = 750;
    b = 1250;

    subplot(1,2,2)
    plot(qgrid(a:b), F1(2,a:b)/F1(2,end), 'linewidth', 2, 'DisplayName','Competitor Distribution: \Delta = \lambda');
    hold on
    plot(qgrid(a:b), F1(end,a:b)/F1(end,end), 'linewidth', 2, 'DisplayName','Competitor Distribution: Marginal');
    hold on
    plot(qgrid(a:b), FNC(a:b)/FNC(end), 'linewidth', 2, 'DisplayName','Non-Competitor Distribution');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Transformed Productivity, q_{hat}", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Conditional Distribution}'}, 'FontWeight', 'normal')
    title("Conditional Productivity Distributions", "FontSize", 14)
    legend('Location', 'southeast')




